clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================


filetosave='Vilcoal_22_04_21_vill_cd_rot'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

village=1;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=3;                         % set to 1 for coalitional deviations, to 2 for self-insurance (SI), to 0 for individual deviations (ID), to 3 for rule of thumb implementation
Tessa1Tobi2=1;                     % sets paths depending on user below
unconditional=0;                   % set this to 1 for unconditional income data, to 0 for income data conditional on time fixed effets / demeaned data period-by-period
one_per_aut=1;                     % set this to 1 if you want deviations to imply one period of autarky
ytaxrate=0;                        % linear tax on income
scatter = 0;		           % set to 1 to just draw a scatter plot
Rov_sb=0;                          % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
                                   % just the utility from consuming average ROV income (first-best)
estrho1=0;			   % set this to 1 if you would like to estimate persistence, default is 0
rhogrid=0;                         % set to 1 if you want to have also a grid for CRRA rho - only to show non-identification in Aurepalle (village=1)
Nphi=1;                            % number of borrowing limits in SI economy - not used currently
server=0;                          % set to 1 to run on server; o wise path for workstations is used
buildup=0;                         % set to 1  to calculate ID model for subset of group sizes, not just whole village
Nincproc=10;                       % number of income processes to look at
if estrho1==0
maxincerr=2/3;                     % maximum income measurment error in multiples of income variance (in levels)
elseif estrho1==1
maxincerr=0;
end
model_var_dy=0;                    % set to 1 if you want to use the model-implied variance of dlog(y) to specify maximum meas error, otherwise var(dlog(y)) from data is used, slightly different as we target var(log(y)), not var(dlog(y())
T=6200; Tinit=201;                 % number of simulation periods in total, T-Tinit+1 equals number of initial periods to discard
maxcerr=0.9;                       % max cons meas error in log levels
if scatter==0
N_cerr=30;                         % number of support points for cons meas error
elseif scatter==1
N_cerr=0;
end
if coaldev==0
    convlength=0;                  % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
else
    convlength=15;
end
if village==1
    vss_high_inddev=34-1;          % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;                
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:7,10,15,vss_high_inddev]; % for buildup==1, this is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1
   
else
    if Tessa1Tobi2==2
        outputpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes/Main Results and Robustness')
        programpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa//JEEA Replication/Programmes/Main Results and Robustness')
    elseif Tessa1Tobi2==1
        outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes/Main Results and Robustness')
        programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
    end
end
        cd(programpath)
if buildup==0
    if rhogrid==0 && scatter == 0
        betavec0=[0.5:0.04:0.78,0.8:0.01:0.98];
elseif scatter ==1
betavec0=[0.83];
    elseif rhogrid==1
        betavec0=[0.83,0.85,0.87,0.89,0.91,0.93,0.95,0.97];
    end
elseif buildup==1

betavec0=[0.86,0.905,0.91,0.94]; %betavec0=[0.7:0.04:0.82,0.84:0.02:0.98];

end
if rhogrid==0
if scatter==1
betavec1=[0.94];
else
 betavec1=[0.84:0.02:0.88,0.89:0.01:0.98];
end
elseif rhogrid==1
    betavec1=[0.84,0.86,0.88,0.90,0.92,0.94,0.96,0.98];
end
betavec3=[0.94:0.005:0.99];
if coaldev<2 || coaldev==3
    
    rhovec=[1];
elseif coaldev==2
    rhovec=ones(1,Nphi);
end

if rhogrid==0
    gapcritset=0.0001;              %convergence criterion for change in values
elseif rhogrid==1
    gapcritset=0.001;
end

rov_average=1;                      % set to 1 to specify the rest of village in terms of average
                                    % consumption
momentclevinlog=1;                  % set this to 1 if you want to calculate the moments for levels of
                                    % consumption based on log consumption
maxitgapset=400;                    % maximum number of iterations
Pun_fac=0.0;                        % consumption penalty in the first period after deviation equal
                                    % to caut*Pun_fac
% grid of time-varying planner weights
lambda=[0.05:0.001:0.95]';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
S_KEEP=nan(1000*3,2,vss_high_inddev);               % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,1000,vss_high_inddev);
PHI_KEEP=nan(m,1000*3,2,vss_high_inddev);
IDIO_DIST_ROV_KEEP=nan(1000*3,3,vss_high_inddev);

% ====================================================
% Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, second: Kanzara, third: Shirapur
if unconditional==0
    datafile='armoments_fe0';
elseif unconditional==1
    datafile='armoments_fe1';
end
if server==0 & Tessa1Tobi2==2
    cd('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Data\Output')
elseif server==0 & Tessa1Tobi2==1
    cd('C:/Users/tbold/Dropbox/Tessa and Tobi/Data/Output')
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);
if model_var_dy==1
    vardy=[0.231,0.125, 0.269];
end
cd(programpath)
if maxincerr>0
    if Nincproc>1
        varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
    else
        varnu_max=zeros(1,3);
    end
    rho_min=covyylag(village)/vary(village);                      % AR(1) coefficient w/o meas error
    rho_max=covyylag(village)/(vary(village)-varnu_max(village)); % AR(1) coefficient w max meas error
    rho1vec=linspace(rho_min,rho_max,Nincproc);
    varnu=max(vary(village)-covyylag(village)./rho1vec,0); % max prevents the first entry to be -e17
    vareps=(vary(village)-varnu).*(1-rho1vec.^2);
else
    rho1_data=covyylag(village)/vary(village);
    vareps_data=vary(village).*(1-rho1_data.^2);
    rho1vec=[-0.6,-0.4,-0.2,0,0.2,0.4,0.6];
    vareps=vary(village).*(1-rho1vec.^2);
    varnu=zeros(1,length(vareps));              
    Nincproc=length(vareps);
end

% ===================================
% Iterate over income processes
% ===================================
    disp('now solving for village,coaldev, unconditional, Nincproc - press any key to continue')
    disp([village coaldev unconditional Nincproc])
    pause
indicator1=zeros(1,size(rhovec,2));
for incproc=1:1%Nincproc
  filetosaveest=['ROT_Momentsest'];

    % specify a markov process for individual income
    [incomevals, Q1] = rouwenhorst(rho1vec(incproc),vareps(incproc)^0.5,2);   
    incomevals=exp(incomevals');
    QQQ=Q1^1000;
    SS=QQQ(1,:)';
    incomevals=incomevals/(SS'*incomevals);
    cumQ=cumsum(Q1,2);
    incomevals=(1-ytaxrate)*incomevals+ytaxrate; % tax as in Krueger and Perri JET 2011
    
    % set minimum and maximum village size
    if coaldev==1                       
        if rhogrid==0
            vss_set=[1:7];% prior solutions have shown that 1-7 is the relevant set. Otherwise use, for example, [1:12,17,22,27,vss_high_coaldev]
        elseif rhogrid==1
            vss_set=[1:5]
        end
        vss_high=vss_set(end);
    elseif coaldev==0
        if buildup==0
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);
        elseif buildup==1
            vss_set=vssbuildupvec_inddev;
            lagind=1;
            vss_high=vss_set(end);
        end
    elseif coaldev==2;
        vss_set=vss_high_inddev;
        vss_high=vss_set(end);
    elseif coaldev==3
        vss_set=[1:10];
        vss_high=10;
    end
    
    
    % ===================================
    % iterate over discount factors
    % ===================================
    if coaldev==0
        betavec=betavec0;
    elseif coaldev==1
        betavec=betavec1;
    elseif coaldev==2
        r=0.04;
if scatter==0
        betamax=1/(1+r)-0.015;
        betavec=linspace(0.8,betamax,20);
elseif scatter ==1
betavec=0.92;
end
        Z=incomevals;
        phinat=min(Z)/(r-1);
        if Nphi>1
            phivec=linspace(phinat,0,Nphi);
        else
            phivec=0;
        end
    elseif  coaldev==3
        betavec=betavec3;
    end
    Momentsest=nan(31,N_cerr+1,size(betavec,2),1,vss_high_inddev);
    for Qbeta=1:size(betavec,2)
        beta=betavec(Qbeta);
        if rhogrid==1
            if coaldev==0
                if beta==0.83
                    rhovec=1.468461538461538
                elseif beta==0.85
                    rhovec=1.316666666666667
                elseif beta==0.87
                    rhovec=1.16
                elseif beta==0.89
                    rhovec=1
                elseif beta==0.91
                    rhovec=0.8340
                elseif beta==0.93
                    rhovec=0.6611
                elseif beta==0.95
                    rhovec=0.4808
                elseif beta==0.97
                    rhovec=0.2912
                end    
            elseif coaldev==1
                if beta == 0.98
                    rhovec=0.366222222222222
                elseif beta==0.84
                    rhovec=2.389
                elseif beta==0.86
                    rhovec=2.11
                elseif beta==0.88
                    rhovec=1.821
                elseif beta==0.90
                    rhovec=1.5515
                elseif beta==0.92
                    rhovec=1.28718
                elseif beta==0.94
                    rhovec=1
                elseif beta==0.96
                    rhovec=0.6935
                else
                end
            end
        end

        % ===================================
        % iterate over CRRA
        % ===================================
        for Qrho=1:size(rhovec,2)
            rho=rhovec(Qrho)  ;
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            
            P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
            waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[];
            WAUT_LAG=nan(1000,1000,length(vss_set));
            STABLE=[];
            for vsscount=1:length(vss_set)       % iterate over village size
                vss=vss_set(vsscount);
                vss_max=max(vss_set);
                if vsscount>1 & coaldev==1
                    lagind=lagind+(vss-(vss_set(vsscount-1)+1));
                end
                disp('now solving for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                % Aiyagari economy
                if coaldev==2
                    sav_problem
                else
                    % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
                    % by vss
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; waut_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov =[]; P_rov =[];
                    S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                    wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
                    Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                    S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];S_rov_aut_all2=[];
                    w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
                    waut=[];
                    
                    % ===================================
                    % build the state space of ROV income
                    % ===================================
                    idio_dist=[];
                    idio_dist_aut=[];
                    idio_dist_aut_all=[];idio_dist_aut_all2=[];
                    idio_dist_sb=[];
                    
                    for k1=0:vss
                        for k2=0:vss
                            
                                II=zeros(1,2);
                                if k1+k2==vss
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals(1,1) + k2*incomevals(2,1)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2];
                                    
                                    % this computes the 'best' distribution of incomes
                                    % of size vss-lagind
                                    if vss>1
                                        II=zeros(1,2);
                                        temp=[k1 k2];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition lagind individuals
                                        % This coalition is one smaller than the
                                        % largest sustainable coalition because it
                                        % will consist of 1 agent and m-1 guys from
                                        % the largest sustainable coalition of size
                                        % m
                                        
                                        for jj=1:lagind
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2]-II;
                                        end
                                        idio_dist_aut=[idio_dist_aut; [k1 k2]-II];
                                        S_rov_aut=[S_rov_aut; ([k1 k2]-II)*incomevals];
                                        
                                        
                                        II=zeros(1,2);
                                        temp=[k1 k2];
                                        % this takes away the lowest incomes to
                                        % derive the distribution of incomes of the
                                        % 'best' coalition
                                        % consisting of m individuals from only the
                                        % rest of village
                                        for jj=1:lagind-1
                                            II(min(find(temp>0)))=II(min(find(temp>0)))+1;
                                            temp=[k1 k2]-II;
                                        end
                                        idio_dist_sb=[idio_dist_sb; [k1 k2]-II];
                                        S_rov_sb=[S_rov_sb; ([k1 k2]-II)*incomevals];
                                    end
                                end
                            
                        end
                    end
                    
                    % resort S_rov in ascending order
                    [S_rov,S_rovind]=sort(S_rov);
                    idio_dist=idio_dist(S_rovind,:);
                    
                    if vss>1
                        idio_dist_aut=idio_dist_aut(S_rovind,:);
                        S_rov_aut=S_rov_aut(S_rovind);
                        idio_dist_sb=idio_dist_sb(S_rovind,:);
                        S_rov_sb=S_rov_sb(S_rovind);
                    end
                    
                    
                    %%%%and create S_rov_aut vectors that drops lagind individuals from
                    %%%%the rov vector. Why lagind? The biggest stable group is
                    %%%%vss+1-lagind, hence, there will be vss+1-lagind-1
                    %%%%invididuals from the rov in it, i.e. vss-lagind, hence
                    %%%%lagind individuals will not be in it. 
                    if vss>1
                        for k=1:length(STABLE)
                            if STABLE(k)==1
                                for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            
                                                if j1+j2==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2) 
                                                    index=index+1;
                                                    temp=[j1 j2 ];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,:,index,k)=[(idio_dist(g,:)-temp)*incomevals];
                                                end
                                            
                                        end
                                    end
                                end
                            end
                            
                        end
                    end
                    
                    % ===================================
                    % simulate aggregate income transitions
                    % ===================================
                    tic
                    T_trans=50000;
                    aggincind=zeros(length(S_rov),length(S_rov));
                    for jjj=1:length(S_rov)
                        jjj;
                        inccount=zeros(T_trans,vss,2);
                        agginc=zeros(T_trans,2);
                        for iii=find(idio_dist(jjj,:)>0)
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ(iii,1)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ(iii,1)<shock & cumQ(iii,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                            end
                            agginc(:,iii)=sum(incomevals(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                        end
                        agginc=sum(agginc,2);
                        for jjj2=1:length(S_rov)
                            qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                            aggincind(jjj,jjj2)=sum(qqq);
                        end
                        P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
                    end
                    disp('Simulation for agg income process takes')
                    toc

                   
                    if vss>1 && coaldev<2
                        S_rov_aut=S_rov_aut/(max(ScaleRov-lagind,1));
                        S_rov_sb=S_rov_sb/(max(ScaleRov-lagind+1,1));
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1
                                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                            end
                        end
                    end
                    S_rov=S_rov/ScaleRov;
                    
                    % ===================================
                    % state space: combinations of indiv and village income
                    % ===================================
                    S=[kron(incomevals,ones(size(S_rov))),kron(ones(size(incomevals)),S_rov)];
                    idio_dist_rov=kron(ones(size(incomevals)),idio_dist);
                    idio_dist_aut=kron(ones(size(incomevals)),idio_dist_aut);
                    S_rov_aut=[kron(ones(size(incomevals)),S_rov_aut)];
                    S_rov_sb=[kron(ones(size(incomevals)),S_rov_sb)];
                    
                    if vss>1 && coaldev<2
                        for g=vss_set(find(vss-1>=vss_set))
                            if STABLE(g)==1 
                                for k=1:size(S_rov_aut_all,3)
                                    S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),S_rov_aut_all(:,:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals)),idio_dist_aut_all(:,:,k,g))];
                                end
                            end
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                    end
                    
                    P=kron(Q1,P_rov);
                    
                    Y=S(:,1)+ScaleRov*S(:,2);			% possible aggregate income outcome
                    state=length(S);                    % number of states
                    agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
                    STATE(vss)=state;
                    
                    
                    % ===================================
                    % Compute set of autarky values
                    % ===================================
                    for k=1:agent
                        caut(:,k)=S(:,k);
                    end
                    waut=zeros(state,agent);
                    waut_avg_lag=zeros(state,agent);
                    waut_sb=zeros(state,agent);
                    % this is just the PDV of utility from consumption of
                    % income caut
                    eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
                    eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
                    if coaldev<2
                                             
                    elseif coaldev==3 % ROT
                        % stationary distribution of aggregate income states
                        SS_S=P^300;
                        SSS_S=SS_S(1,:);
                        SSS_Ssave(vss+1,1:length(SSS_S))=SSS_S;
                        sizeSSS_Ssave(vss+1)=length(SSS_S);
                        SOI=[0.5 0.6 0.7 0.8 0.9 1]  ; %%strenght of insurance, i.e. what share of first-best transfer is made
                        QUASI=[0 0.025 0.05]; %%quasi-lending, i.e. what share of previous transfers is repaid in periods where all have same income realization
                        % value from consuming average income
                        if vss==1  %here we are definining the autarky payoffs. 
                            sizeSSS_Ssave(1)=2;
                            SSS_Ssave(1,1:2)=SS;
                            
                            V_ROT(1,1:2,1)=(inv(diag(ones(length(incomevals),1))-beta*Q1)*utility(rho,incomevals));
                            V_ROT(1,1:2,2)=(inv(diag(ones(length(incomevals),1))-beta*Q1)*utility(rho,incomevals));
                            P_LAG(1:2,1:2,1)=Q1;
                            Y_lag(1,1:2,1)=incomevals';
                            
                            for kkk=1:length(SOI)
                            for ggg=1:length(QUASI)
                            V_ROT(1,1:2,1:2,kkk,ggg)=V_ROT(1,1:2,1:2,1,1);
                            end
                            end
                            
                        end
                        
                        
                        
                                         
                        % vvss are the number of high income individuals
                        Y_lag(vss+1,1:length(Y))=Y;
                        P_LAG(1:length(S),1:length(S),vss+1)=P;
                        S_rov_KEEP(1:length(S)/2,vss)=S_rov;
                        S_KEEP(1:length(S),:,vss)=S;
                        IDIO_DIST_ROV_KEEP(1:length(S),1:size(idio_dist_rov,2),vss)=idio_dist_rov;
                        %%%%and now we do the same thing with some history
                        %%%%dependence in periods where everyone has the
                        %%%%same income               
                        
                        %this calculates the per period utility for the
                        %individual when they are in a group of size vss+1    
                        for kkk=1:length(SOI)
                        for ggg=1:length(QUASI)
                        quasi=QUASI(ggg);
                        soi=SOI(kkk);
                        vv=zeros(state,state);
                        for tt=1:5000
                        for j1=1:state
                        for j2=1:state
                        if (j1==1 | j1==state) && (j2~=1 & j2~=state) && S(j2,1)==max(incomevals)
                        transfer=-soi*(min(incomevals)-Y(j2)/(vss+1)); %now this is the transfer the LI individual received
                        R(j1,j2)=utility(rho,S(j1,1)+idio_dist_rov(j2,1)/(idio_dist_rov(j2,2)+1)*quasi*transfer);
                        elseif (j1==1 | j1==state) && (j2~=1 & j2~=state) && S(j2,1)==min(incomevals)
                        transfer=-soi*(min(incomevals)-Y(j2)/(vss+1)); %now this is the transfer the LI individual received 
                        R(j1,j2)=utility(rho,S(j1,1)-quasi*transfer);
                        elseif (j1~=1 & j1~=state) && S(j1,1)==min(incomevals)
                        R(j1,j2)=utility(rho,S(j1,1)-soi*(S(j1,1)-Y(j1)/(vss+1)));
                        elseif (j1~=1 & j1~=state) && S(j1,1)==max(incomevals)
                        R(j1,j2)=utility(rho,S(j1,1)-soi*(S(j1,1)-Y(j1)/(vss+1)));
                        elseif (j1==1 | j1==state) && (j2==1 | j2==state)
                        R(j1,j2)=utility(rho,S(j1,1));   
                        end
                            
                        vv(j1,j2)=R(j1,j2)+beta*P(j1,:)*[vv(:,j1)];
                        
                        vv(2,3);
                        end
                        end
                        end
                        V_ROT(vss+1,(1:size(vv,1)),1:size(vv,2),kkk,ggg)=vv; 
                        %%%%%%%%participation constraints
                        end
                        end
                        if vss>1
                            
                            for nnn=1:vss
                                for kkkdev=1:length(SOI)
                                    for gggdev=1:length(QUASI)
                                        if nnn>1
                                        if STABLE_ROT(nnn-1,kkkdev,gggdev)==nnn-1
                                        VDEV(kkkdev,gggdev,nnn)=P_LAG(sizeSSS_Ssave(nnn),1:sizeSSS_Ssave(nnn),nnn)*V_ROT(nnn,1:sizeSSS_Ssave(nnn),sizeSSS_Ssave(nnn),kkkdev,gggdev)';
                                        else
                                        VDEV(kkkdev,gggdev,nnn)=-9999;   
                                        end
                                        elseif nnn==1
                                        VDEV(kkkdev,gggdev,nnn)=P_LAG(sizeSSS_Ssave(nnn),1:sizeSSS_Ssave(nnn),nnn)*V_ROT(nnn,1:sizeSSS_Ssave(nnn),sizeSSS_Ssave(nnn),kkkdev,gggdev)';
                                        end
                                    end
                                end
                            end
                            for nnn=1:vss
                            VDEV_RS=reshape(VDEV(:,:,1:nnn),size(VDEV,1)*size(VDEV,2)*nnn,1);
                            [CC,II]=max(VDEV_RS);
                            [i1, i2, i3]=ind2sub(size(VDEV(:,:,1:nnn)),II)    ;
                            kkkmax(nnn)=i1;
                            gggmax(nnn)=i2;
                            nnnmax(nnn)=i3;
                           
                            end
                            clear VDEV VDEV_RS 
                        else 
                            kkkmax(1)=1;
                            gggmax(1)=1;
                            nnnmax(1)=1;
                        
                        end
                        
                        for kkk=1:length(SOI)
                        soi=SOI(kkk);
                        for ggg=1:length(QUASI)
                        for vvss=1:vss
                        
                        ind1=find(abs(Y-(vvss)*max(incomevals)-(vss+1-vvss)*min(incomevals))<0.00000001 & S(:,1)==max(incomevals)); % this is the indicator of the state today with vvss high and vss+1-vvss low income individuals - can be several
                            % this is the utility loss from making transfers today when
                            % vvss out of vss individuals have high income and vss-vvss
                            % have minimum income
                            deltav(vss,vvss)=(utility(rho,max(incomevals))-utility(rho,(max(incomevals)-soi*(max(incomevals) - (vvss*max(incomevals)+(vss+1-vvss)*min(incomevals))/(vss+1)))));
                            %v(vss,vss)=(1-delta)* (utility(gamma,max(incomevals))+delta*SSSsave(vvss,1:sizeSSSSsave)*P(1,:)*V(vvss);
                            % the difference in continuatoin values from having full
                            % insurance with vss vs vvss individuals
                                                        
                                               
                            deltaV(vss,vvss)=beta*(P(ind1(1),:)*V_ROT(vss+1,1:sizeSSS_Ssave(vss+1),ind1(1),kkk,ggg)'-P_LAG(sizeSSS_Ssave(nnnmax(vvss)),1:sizeSSS_Ssave(nnnmax(vvss)),nnnmax(vvss))*V_ROT(nnnmax(vvss),1:sizeSSS_Ssave(nnnmax(vvss)),sizeSSS_Ssave(nnnmax(vvss)),kkkmax(vvss),gggmax(vvss))');
                            
                            % if this is negative for vvss, and the group with vvss is stable, then full insurance with vss individuals
                            % is sustainable against a deviation with vvss individuals
                            Delta(vss,vvss)=deltav(vss,vvss)-deltaV(vss,vvss);
                            
                        end
                        stable=find(Delta(vss,:)<0)
                        if length(stable)==vss
                        kkk
                        STABLE_ROT(vss,kkk,ggg)=vss;
                        else
                        STABLE_ROT(vss,kkk,ggg)=0;
                        end
                        end
                        end
                        
                        end
                        end
                        
                        
                        clear yind  y yy dy aggyy aggy 
                        
                        
                        end
                
                        
                   
                
                
            
            
            
            %==========================================================
            % Simulate the economy
            % =========================================================
            % This lists all possible combinations of individual income
            % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
            % with every individual income value
for k_soi=1:length(SOI)
    for g_quasi=1:length(QUASI)
            if coaldev==1                       
                stable=(find(STABLE==1));
                vss_set_sim=stable;
            elseif coaldev==0                             
                vss_set_sim=vss_set;
            elseif coaldev==2
                vss_set_sim=vss_set;
            elseif coaldev==3
                stable=max(find(STABLE_ROT(:,k_soi,g_quasi)~=0)); %%%only do the simulation for the maximum stable group size
                vss_set_sim=stable;
            end

            for vss=[vss_set_sim]
                              
                Nvill=ceil((vss_high_inddev+1)/(vss+1))
                var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
                varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
                vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
                vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
                R2_dcdy=zeros(1,Nvill*(vss+1));betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
                beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
                beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
                
                %=============================================
                % Start simulation
                % ============================================
                % number of insurance groups equals village size
                % vss_high+1 divided by size of ins group vss+1
                
                Npanel=Nvill*(vss+1);
                csim0=nan(Nvill,vss+1,T);
                income0=nan(Nvill,vss+1,T);
                income0ind=nan(Nvill,vss+1,T);
                phisim0=nan(Nvill,vss+1,T);
                Yvill=nan(1,T);
                Yagg=nan(Nvill,T);
                Ssim=nan(Nvill,vss+1,T);
                gammarfun=nan(Nvill,vss+1,T);
                phinew=nan(Nvill,vss+1,T);
                indSjj=nan(T,vss+1);
                disp('simulating') 
                disp('Reminder: now simulating for village,coaldev,incproc')
                disp([village coaldev incproc])
                disp('now simulating for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                tic
                %initial state
                t=1;
                % initial values of income and consumption (set =
                % income)
                rng(3,'twister');
                shockinittemp=randi(2,Nvill*(vss+1),1);
                shockinit=reshape(shockinittemp,Nvill,vss+1);
                income0ind(:,:,1)=shockinit;
                income0(:,:,1)=incomevals(income0ind(:,:,1));
                csim0(:,:,1)=income0(:,:,1);
                Yagg(:,t)=sum(income0(:,:,t),2);
                IDIO_DIST_sim=[];
                
                while t<T
                    t=t+1;
                    %%%%simulate next period income taking into account
                    %%%% persistence across time
                     rng(t,'twister');
                    shocktemp=rand(Nvill*(vss+1),1);
                    shock=reshape(shocktemp,Nvill,vss+1);
                    for k=1:vss+1
                        rr=(find(shock(:,k)<=cumQ(income0ind(:,k,t-1),1)));
                        income0ind(rr,k,t)=1;
                        rr=(find(cumQ(income0ind(:,k,t-1),1)<shock(:,k) & cumQ(income0ind(:,k,t-1),2)>=shock(:,k)));
                        income0ind(rr,k,t)=2;
                        income0(:,k,t)=incomevals(income0ind(:,k,t))  ;
                    end
                    Yagg(:,t)=sum(income0(:,:,t),2);
                end
                
                if coaldev<2
                 
                elseif coaldev==2
                
                elseif coaldev==3
                    
                   state=STATE(vss);
                   t=0
                   while t<T
                       t=t+1;
                   for ivill=1:Nvill
                      %%%find the correct states
                    
                                           
                      for k=1:vss+1
                      %this gives the no of high income inds except
                      %individual under consideration
                      rrhigh_j1=sum(income0ind(ivill,:,t)==2);
                      rrlow_j1=sum(income0ind(ivill,:,t)==1);
                      if t>1
                      rrhigh_j2=sum(income0ind(ivill,:,t-1)==2);
                      rrlow_j2=sum(income0ind(ivill,:,t-1)==1);
                      end
                      if income0ind(ivill,k,t)==2
                      rrhigh_j1=rrhigh_j1-1;
                      else
                      rrlow_j1=rrlow_j1-1;
                      end
                      if t>1 && income0ind(ivill,k,t-1)==2
                      rrhigh_j2=rrhigh_j2-1;
                      elseif t>1 && income0ind(ivill,k,t-1)==1
                      rrlow_j2=rrlow_j2-1;
                      end
                    j1=find(abs(S_KEEP(1:state,1,vss)-income0(ivill,k,t))<0.0001 ...
                               & rrlow_j1==IDIO_DIST_ROV_KEEP(1:state,1,vss) & rrhigh_j1==IDIO_DIST_ROV_KEEP(1:state,2,vss));
                    j2=1;
                    if t>1
                    j2=find(abs(S_KEEP(1:state,1,vss)-income0(ivill,k,t-1))<0.0001 ...
                               & rrlow_j2==IDIO_DIST_ROV_KEEP(1:state,1,vss) & rrhigh_j2==IDIO_DIST_ROV_KEEP(1:state,2,vss));
                    end
                    soi=SOI(k_soi);
                    quasi=QUASI(g_quasi);
                    if (j1==1 | j1==state) && (j2~=1 & j2~=state) && S_KEEP(j2,1,vss)==max(incomevals) && t>1
                        transfer=-soi*(min(incomevals)-Yagg(ivill,t-1)/(vss+1));
                        csim0(ivill,k,t)=S_KEEP(j1,1,vss)+rrlow_j2/(rrhigh_j2+1)*quasi*transfer;
                        elseif (j1==1 | j1==state) && (j2~=1 & j2~=state) && S_KEEP(j2,1,vss)==min(incomevals) && t>1
                        transfer=-soi*(min(incomevals)-Yagg(ivill,t-1)/(vss+1));
                        csim0(ivill,k,t)=S_KEEP(j1,1,vss)-quasi*transfer;
                        else
                        csim0(ivill,k,t)=S_KEEP(j1,1,vss)-soi*(S_KEEP(j1,1,vss)-Yagg(ivill,t)/(vss+1));
                        end
                    
                       end
                   end
                   end
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                end
                disp('simulating takes')
                tsim(Qbeta,Qrho)=toc;            
                disp(tsim(Qbeta,Qrho))
                
                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                 rng(3,'twister'); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
                 rng(4,'twister'); incerr=randn(size(csim0));   Momentsest1=nan(30,N_cerr); ind=[];
                % loop over size of cons measurement error
                for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                for ic=1:N_cerr+1
                    if N_cerr>0
                        csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    else
                        csim=exp(log(csim0));
                    end
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);
                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    end
                    aggincome=sum(income,1);
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); 
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                    
                    % calculate moments for each individiual
                    for ind=1:size(csim,1)
                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                        
                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease!!
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%%conditional on village income
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
                        
                        %%%conditional on group income
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                      
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%% Relative variances
                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
                        
                        %%% THE BETA COEFFICIENTS
                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc(ind,:)',[ones(1,size(yy,2))',aggyy(1,:)',yy(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt(ind,vss)=betadcdytemp(2);
                        RR2_dcdy(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
                        % for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(:,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        
                        %2) conditional on village income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);

                        %3) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    R2_dcdy=mean(RR2_dcdy(1:size(csim,1),vss));
                    
                    % =====================
                    % Matrix of moments
                    % =====================                    
                    % col 1-5 are parameters
                    Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                    mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                    mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
                    Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
                    Momentsest1(31,ic)=mean(RR2_dcdy(1:size(csim,1),vss));
                    %%% and some other model-dependent measures
                    if coaldev<2
                        Momentsest1(29,ic)=CRITGAP(vss);
                    elseif coaldev==2
                        Momentsest1(29,ic)=phi;
                    elseif coaldev==3
                        Momentsest1(29,ic)=nan;
                    end
                end
                
                Momentsest(:,:,Qbeta,Qrho,vss,k_soi,g_quasi)=[Momentsest1];
                
                toc
            end
        end
    end

                cd(outputpath)
                save(filetosaveest,'Momentsest');
                cd(programpath)
        end
    end
    
    cd(outputpath)
   
        eval(['save ', filetosave '.mat'])
    
    cd(programpath)
end